source('00_util_scripts/mod_bplot.R')
library(MSstats)

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

cytof_mono <-
read_csv('mission/SLE_TRPM2_MfMo/data/combat22sepsis/cytof/CBD-KEY-CYTOF-MYELOID/monocyte_expression.csv')

cytof_mono$cluster_id |> unique()

cytof_mono |>
  filter(cluster_id == 'ki67') |>
  ggplot(aes(condition, Freq)) +
  geom_boxplot()

sg_time <- read_csv('searchgui.csv')

sg_shake <- read_tsv('searchgui_shake.csv')

sg_shake |>
  pivot_longer(-1) |>
  mutate(Engine = fct_reorder(Engine, value, max)) |>
  ggplot(aes(value, Engine, fill = name)) +
  geom_col() +
  facet_wrap(~str_to_title(name), scales = 'free_x') +
  labs(x = 'Protein count', title = 'PeptideShaker protein recovery') +
  theme_pubr(legend = 'none')

sg_shake |> left_join(sg_time) |>
  filter(total > 0) |>
  pivot_longer(2:4) |>
  mutate(effiency = value/Time_Seconds) |>
  ggplot(aes(effiency, Engine)) +
  geom_col() +
  facet_wrap(~name, scales = 'free_x')

# filter for sepsis/HV timsTOF raw data -------------
timstof_ftp <- read_tsv('mission/SLE_TRPM2_MfMo/data/combat22sepsis/timsTOF_ftp.tsv')

timstof_ftp |>
  filter(str_detect(NAME, 'sepsis|HV'), TYPE == 'RAW') |>
  select(URI) |>
  write_csv('combat_timsTOF_dzip.uri', col_names = F)

timstof_md5 <- timstof_ftp |>
  filter(NAME == 'checksum.txt') |>
  pull(URI) |>
  read_tsv(skip = 1, col_names = F)

timstof_md5 |>
  mutate(X1 = str_remove(X1, '.+\\\\')) |>
  write_tsv('combat_timsTOF_md5sum.txt',col_names = F)

# fragpipe 2 msstats -----------
fp_pepsis <-
  read_tsv('mission/SLE_TRPM2_MfMo/data/combat22sepsis/combined_protein.tsv')

fp_msstat <-
  data.table::fread('~/append-ssd/data_lfs/msstats.csv.gz')

ms_pepsis <- FragPipetoMSstatsFormat(fp_msstat)

# take 3min
sum_pepsis <- ms_pepsis |>
  dataProcess(remove50missing = T, numberOfCores = 4)

model <- groupComparison("pairwise", sum_pepsis)

sepsis_dep <- model$ComparisonResult |>
  as_tibble() |>
  separate_wider_delim(Protein, delim = '|', names = c(NA, 'UNIPROT', NA),
                       cols_remove = F)

model$ModelQC |>
  as_tibble() |>
  slice_max(ABUNDANCE)

sepsis_dep <- sepsis_dep$UNIPROT |>
  clusterProfiler::bitr(fromType = 'UNIPROT', toType = 'SYMBOL',
                        OrgDb = 'org.Hs.eg.db') |>
  left_join(sepsis_dep)

sepsis_dep |>
  mutate(avg_log2FC = -log2FC,
         p_val_adj = ifelse(adj.pvalue == 0, 1e-10, adj.pvalue), gene = SYMBOL) |>
  filter(adj.pvalue != 0) |>
  plot_bill_volc(group1 = 'Sepsis', group2 = 'HC') +
  labs(title = 'Differential enriched plasma protein in sepsis patients',
       subtitle = 'PXD023175 (62 Sepsis vs 30 HC)')

sepsis_dep |>
  write_source_csv('combat.sepsis.plasma.protein.logfc')

outlier <- sepsis_dep |>
  filter(adj.pvalue == 0) |>
  distinct(Protein, .keep_all = T)

outlier_logint <- sum_pepsis |>
  quantification(format = 'long', use_log_file = F) |>
  as_tibble() |>
  right_join(outlier)

outlier_logint |>
  mutate(LogIntensity = ifelse(is.na(LogIntensity), 1, LogIntensity),
         group = str_extract(Group_Subject, 'HC|sepsis')) |>
  ggplot(aes(group, LogIntensity)) +
  geom_boxplot() +
  geom_jitter() +
  facet_wrap(~Protein, scales = 'free_y')

outlier_logint |>
  mutate(LogIntensity = ifelse(is.na(LogIntensity), 1, LogIntensity),
         group = str_extract(Group_Subject, 'HC|sepsis'),
         SYMBOL = fct_reorder(SYMBOL, log2FC)) |>
  ggplot(aes(Group_Subject, SYMBOL, fill = LogIntensity)) +
  geom_raster() +
  theme_pubr(legend = 'right') +
  scale_fill_viridis_c() +
  theme(axis.text.x = element_blank()) +
  labs(x = '', y = 'Protein', title = 'Sepsis vs HC plasma protein')

g1 <- last_plot()  

g2 <- outlier_logint |>
  mutate(group = str_extract(Group_Subject, 'HC|sepsis'),
         protein = str_remove_all(Protein, '.+\\||_HUMAN') |> fct_reorder(log2FC)) |>
  ggplot(aes(Group_Subject, '', fill = group)) +
  geom_raster() +
  theme_void() +
  scale_fill_hue(direction = -1) +
  theme(axis.text.x = element_blank())
  
g1 / g2 + patchwork::plot_layout(heights = c(15,1), guides = 'collect')
